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^Abstract 

^ombining fast MR acquisition sequences and high resolution imaging is a major issue in 
Cjlynamic imaging. Reducing the acquisition time can be achieved by using non-Cartesian 
and sparse acquisitions. The reconstruction of MR images from these measurements is 
^^enerally carried out using gridding that interpolates the missing data to obtain a dense 
^Cartesian /c-space filling. The MR image is then reconstructed using a conventional Fast 
l^'ourier Transform (FFT). The estimation of the missing data unavoidably introduces 
C^irtifacts in the image that remain difficult to quantify. 

m ' 

• A general reconstruction method is proposed to take into account these limitations. It can 
i^ye applied to any sampling trajectory in fc-space, Cartesian or not, and specifically takes 
Ointo account the exact location of the measured data, without making any interpolation of 
( ^h'e missing data in fc-space. Information about the expected characteristics of the imaged 
Object is introduced to preserve the spatial resolution and improve the signal to noise 
^■atio in a regularization framework. The reconstructed image is obtained by minimizing a 
^non-quadratic convex objective function. An original rewriting of this criterion is shown 
" to strongly improve the reconstruction efficiency. Results on simulated data and on a real 
spiral acquisition are presented and discussed. 

Key words: Fast MRI, Fourier synthesis, inverse problems, regularization, 
edge-preservation. 



1 Introduction 



In Magnetic Resonance Imaging (MRI) the acquired 
data are samples of the Fourier transform of the im- 
aged object [1]. Acquisition is often discussed in terms 
of location in A;-space and most conventional methods 
collect data on a regular Cartesian grid. This allows 
for a straightforward characterization of aliasing and 
Gibbs artifacts, and permits direct image reconstruc- 



tion by means of 2D-Fast Fourier Transform (FFT) al- 
gorithms. Other acquisition sequences, such as spiral 
[2], PROPELLER [3], projection reconstruction, i.e. 
radial [4], rosette [5], collect data on a non-Cartesian 
grid. They possess many desirable properties, includ- 
ing reduction of the acquisition time and of various 
motion artifacts. The gridding procedure associated to 
an FFT is the most common method for Cartesian im- 
age reconstruction from such irregular A;-space acqui- 
sitions. 



Re-gridding data from non-Cartesian locations to a 
Cartesian grid has been addressed by many authors. 
O' Sullivan [6] introduced a convolution-interpolation 
technique in computerized tomography (CT) which 
can be applied to magnetic resonance imaging [2]. He 
suggested not to use a direct reconstruction, but to per- 
form a convolution-interpolation of the data sampled 
on a polar pattern onto a Cartesian /c-space. The final 
image was obtained by FFT. The stressed advantage 
of this technique was the reduction of computational 
complexity compared to the filtered back-projection 
technique. Moreover, it can be applied to any arbitrary 
trajectory in /c-space. 

More generally, the reconstruction process is four 
steps: 

(1) data weighting for nonuniform sampling com- 
pensation, 

(2) re-sampling onto a Cartesian grid, using a given 
kernel, 

(3) computation of the FFT, 

(4) correction for the kernel apodization. 

Jackson et al. [7] precisely discussed criteria to choose 
an appropriate convolution kernel. This is necessary 
for accurate interpolation and also for minimization 
of reconstruction errors due to uneven weighting of 
/c-space. Several authors have suggested methods for 
calculating this sampling density. Numerical solutions 
have been proposed that iteratively calculate the com- 
pensation weights [3]. But, for arbitrary trajectories, 
the weighting function is not known analytically and 
must somehow be extracted from the sampling func- 
tion itself. A possible solution is to use the area of the 
Voronoi cell around each sample [8]. 

The gridding method is computationally efficient. 
However, convolution-interpolation methods unavoid- 
ably introduce artifacts in the reconstructed images 
[8]. Indeed, for a given kernel the convolution modi- 
fies data in A;- space and it is difficult to know the exact 
effect of gridding in the image domain. Moreover, this 
method tends to correlate the noise in the measured 
samples and lacks solid analysis and design tools to 
quantify or minimize the reconstruction errors. 

The principle of regularized reconstruction has been 
described by several authors for parallel imaging: [9], 
[10] and more recently [11] proposed the use of a 



general reconstruction method for sensitivity encod- 
ing (SENSE) [12] which has been applied with a 
quadratic regularization term and a Cartesian acquisi- 
tion scheme. In this paper, we extend this work by: 1) 
giving a more general formulation of the reconstruc- 
tion term for Non Cartesian trajectories, 2) specifically 
using the exact non-uniform locations of the acquired 
data in A> space, without the need for gridding the data 
to a uniform Cartesian grid and, 3) incorporate a non- 
quadratic convex regularization term in order to main- 
tain edge sharpness compared to a purely quadratic 
term. The regularization term represents the prior in- 
formation about the imaged object that improves the 
signal to noise ratio (SNR) of the reconstructed image 
as well as the spatial resolution. 

In section 2, we recall the basis of MRI signal ac- 
quisition and the modelling of the MR acquisition 
process. Then we address the image reconstruction 
methods for different acquisition schemes and develop 
the proposed method, in section 3. The reconstruc- 
tion is based on the iterative optimization of a Dis- 
crete Fourier Transform (DFT) regularized criterion. 
Rewriting this criterion allows to reduce the complex- 
ity of the computation and to decrease the reconstruc- 
tion time. Finally, section 4 compares the proposed 
method and the gridding reconstruction for simulated 
and real sparse data acquired along interleaved spiral 
trajectories. 



2 Direct model 

MRI theory [1] indicates that the acquired signal s is 
related to the imaged object / through: 

*(*(*)) = JJ D f(r)e^ tr dr, (1) 

in a 2D context. D is the field of view, i.e., the extent 
of the imaged object, r is the spatial vector and k(t) = 
[k x (t), kyft)] 1 ("t" denotes a transpose) is the A;-space 
trajectory. Thus, the received signal can be thought as 
the Fourier transform of the object, along a trajectory 
k(t) determined by the magnetic gradient field G(t) = 

k(t) = 7 r Git') at'. 

Jo 



The modulus of f(r) is proportional to the spin density 
function and the phase factor is influenced by spin 
motions and magnetic field inhomogeneities. 

Remark 1 — Eq. (1) presents a model for an ideal 
signal. Actual signals also include terms for the relax- 
ation of the magnetic moments which will cause the 
signal amplitude to decrease, as well as a term for in- 
homogeneity within the image. By the way, they could 
be easily incorporated in (1), but for our purposes 
here we will ignore these effects. 

Practically, the acquired signal is not a continuous 
function of time but made of a finite number of sam- 
ples. This introduces the discretization of the data, and 
the measured data set writes s = [s , s±, . . . , sl-i]* £ 
C L , i.e., consists of L data sampled along the discrete 
trajectory [k , k 1} . . . , k L -i], where ki = [k l x , k$. For 
a single sample, Eq. (1) then reads: 



si 



' dr . 



Generally the object / is not reconstructed as a con- 
tinuous function of the spatial variables r but is also 
discretized for practical considerations: to use image 
visualization and also to perform fast reconstruction 
techniques by means of FFT This introduces a dis- 
cretization of the unknown object and a common 
choice is a Cartesian grid of size N xN. We note f rhm 
the unknown discretized object evaluated at locations 
fnm = [n, to]* with n, m = 0, 1, . . . , N — 1. 

The discrete model is then given by an approximation 
of the integral of Eq. (1): 



si 



1 



N-l 

v 

iv n,m=0 



fn,r, 



j2n(k l x m/F x +h l y n/F y ) 



where F = [F x , Fyf is the spatial sampling frequency 
of the object. To comply with the Shannon sampling 
frequency, F must be chosen such as F x > 2/D x and 
F y > 2/Dy, where D x and D y are the dimensions of 
the field of view. For sake of simplicity we assume 
here that F = [1,1]* and the spatial frequencies k l x 
and k l y are normalized and lie in [—0.5, +0.5]. 

In practice the acquired samples are corrupted by a 
complex valued noise, denoted b = [fe , • • • , 61,-1]* £ 
C L , which can be assumed to be additive white and 
Gaussian [13]. 



We can then write, for one datum, the final discretized 
model as: 



si = 



1 

N 



N-l 



£ Um^^^+h (2) 
n,m=0 

for / = 0, . . . , L — 1 or, more simply as 

si = hj + bi , 

with / being a column vector, collecting the f n>m re- 
arranged column by column in one vector, and h t a 
row vector 



^ _ J_ ri27rfc*T-oo e i2irfc*r i . Jvi, ; ■•• , - ,] 



N 



, . . . , e 



The whole data vector then writes: 

s = Hf + b, 
where H is the inverse Fourier matrix: 



(3) 



H = 



\h L -i) 

depending on the acquisition locations. 

Eq. (3) is a linear model with additive Gaussian noise. 
It has been extensively studied in literature [14]. The 
aim of the reconstruction process is to compute an es- 
timate / of the unknown object / from the discrete, 
incomplete and noisy /c-space samples s. The prob- 
lem is referred to as a Fourier synthesis problem and 
consists of inversion of the model (3). 



3 Model inversion 



A usual inversion method relies on a Least Squares 
(LS) criterion, based on Eq. (3): 

JLsCf) = ||* - Hf\\ 2 = J2 \ Sl - hj\ 2 . (4) 

1=0 

The reconstructed image is the minimizer of J LS : 

/ LS = argmin J LS (f) , 
/ 



and minimizes the quadratic error between the mea- 
sured data and the estimated ones generated by the 
direct model (3). The solution writes: 

/ LS = (^)- 1 frt a> 

if H^H is invertible, property that depends on the ac- 
quisition scheme. 

3.1 Cartesian and complete acquisitions 

In Complete Cartesian ( CC) acquisitions H is the N x 
N inverse Fourier transform matrix, evaluated on an 
uniform grid. We then have WH = I and the LS 
solution simplifies to 

f = H*8. (5) 

It is efficiently computed by the FFT of the raw data 
and the compromise between acquisition time and im- 
age characteristics depends only on the acquisition 
scheme. 

This inversion method directly holds as long as a com- 
plete Cartesian /c-space is available as for the conven- 
tional line by line acquisitions where one line is ac- 
quired for each successive radio-frequency (rf) exci- 
tation. It holds also for multi-shot acquisitions when 
more than a single A;-space line is acquired for each rf 
excitation. It can finally be applied to EPI sequences 
when only one excitation is used to sample the whole 
/c-space domain. 

The method remains convenient for time segmented 
acquisitions that update only partially A>space, such as 
keyhole, BRISK or TRICKS techniques [15,16,17,18] 
provided that a convenient filing of /c-space data has 
been made previously. 

3.2 Incomplete and non Cartesian acquisitions 

Other acquisition schemes have been proposed in or- 
der to reduce acquisition time. They can be divided in 
two groups: Incomplete Cartesian (IC) ones and Non 
Cartesian (NC) ones. 

IC - Partial Cartesian filling of a k- space such as the 



widely used "half Fourier" method [19] or vari- 
able density phase encoding technique [20] allow 
to reduce the number of acquired data and thus 
the acquisition time. In this case, H is a partial 
matrix and can still be computed with the FFT. 
NC - Non Cartesian /c-space filling (interleaved spi- 
rals, PROPELLER sequence, radial, concentric 
circles, rosettes...) conjugate a variable, non- 
uniform density encoding with specific gradient 
sequences with the same objective of acquisi- 
tion time reduction. These acquisition schemes 
often require a small number of rf pulses, take 
advantage of the available gradient strength and 
rising time, reduce motion artifacts and lessen 
sensitivity to off-resonances and field inhomo- 
geneities [2]. 

From a mathematical stand point, the main difficulty 
of the Non Cartesian acquisition schemes is that (5) 
cannot be computed using the FFT algorithm, since 
the samples are no longer on a uniform grid. Current 
strategies force the re-use of FFT reconstruction (5) 
by means of data pre-processing. 

IC - The missing data are completed beforehand using 
Fourier symmetry properties of the A:- space [19] 
(see also the Margosian reconstruction [21]), or 
a zero-padding extrapolation. Conventional zero 
padding used to construct a square image from 
a rectangular acquisition matrix also belongs to 
this category. 

NC - The acquired data are interpolated and re- 
sampled by means of a gridding method. 

Thus a complete Cartesian A;-space is pre-computed 
from the acquired data and the final image is obtained 
by FFT. The wide availability of high-speed FFT rou- 
tines and processors have made the method by far the 
most popular. But, such methods do not rely on the 
physical model (3) nor on the true acquired data: they 
introduce interpolated data resulting in inaccuracies 
in the reconstructed images. On the contrary, the pro- 
posed method accounts for exact locations of the data 
in /c-space. The methodology is applicable for both IC 
and NC acquisition scheme and we concentrate on the 
NC case i.e. the non-uniform DFT model. 

Other strategies rely on true DFT and LS framework. 
The main problem here is that H^H is not invert- 
ible: the unknown image pixels usually outnumber 



the acquired data and the problem is indeterminate, 
i.e., J^ s does not have a unique minimizer. From ba- 
sic inverse problem theory, several regularization ap- 
proaches have been proposed. Among the earliest are 
the Truncated Singular Value Decomposition (TSVD) 
and the Minimum Norm Least Squares (MNLS). They 
properly regularize the problem, alleviate the indeter- 
minacy and define a solution to (3). The TSVD and the 
MNLS approaches have been proposed in MRI by [20] 
for IC acquisition and by [22,23] for NC acquisitions, 
respectively. Practically they both can be extended for 
IC and NC acquisitions and behave similarly. 

In any case (TSVD, MNLS, gridding, zero-padding), 
it is difficult to control the information accounted for, 
in order to regularize the problem. Moreover they 
cannot incorporate more specific information such as 
pixel correlation, and edge enhancement. The pro- 
posed method, described below, accounts for known 
common information about the expected images and 
exact locations of the data in the /c-space. 



3.3 Regularized Method 

The proposed method relies on Regularized Least 
Squares (RLS) criterion: 

It is based on the LS term and a prior one 1Z, that 
only depends upon the object /. The proposed solution 
writes: 

/ Reg = argmin J Ree (f) ■ 
f 

The choice of 1Z depends on the information to be in- 
troduced. In MR, there are a great variety of image 
kinds, but at least two common characteristics are ob- 
served. 

(1) The structure have usually smooth variations and 
a good contrast compared to the surrounding or- 
gans, more particularly when contrast agents are 
used. These regions are separated by sharp tran- 
sitions representing the edges. 

(2) The regions outside the imaged object i.e. the 
background is a region where / is expected to 
be zero. 



The proposed regularization term accounts for these 
information and takes the following form: 

iz(f) = \ 1 n 1 (f) + \ n (f). 

A x and A are the regularization parameters (hyper- 
parameters) that balance the trade-off between the fit 
to the data and the prior. One can clearly see that 
A x = A = gives the LS criterion, and no informa- 
tion about the object is accounted for. On the contrary, 
when Ax, A — > oo the solution is only based on the a 
priori information. 

The first term fix (/) is an edge-preserving smoothness 
term based on the first order pixel differences in the 
two spatial directions: 

— X! W (fn+l,m — fn,m) 
n,m 

,m+l fn,m) i 

n,m 

and the second one fio(/) introduces the penalization 
for the image background: 

fi o(/) = $Z^«o(/n,m) . 

n,m 

The penalization functions (p a parametrized by the co- 
efficient a (discussed below) determine the character- 
istics of the reconstruction and has been addressed by 
many authors [24,25,26,27,28]. 

Interesting edge-preserving functions are those with 
a flat asymptotic behaviour towards infinity, such as 
the Blake and Zisserman function [27] or Geman and 
McClure [28]. However these functions are not con- 
vex and the resulting regularized criterion may present 
numerous local minima. Its optimization therefore re- 
quires complex and time-consuming techniques. On 
the contrary, the quadratic function proposed by Hunt 
[25]: (p(x) = x 2 is best suited to fast optimization 
algorithms. Nevertheless, it tends to introduce strong 
penalizations for large transitions (see Fig. 1), which 
may over-smooth discontinuities. An interesting trade- 
off can be achieved by using a combination between 
a quadratic function (L 2 ) to smooth small pixel dif- 
ferences and a linear function (Li) for large pixel 
differences beyond a defined threshold a. The latter 
part produces a lower penalization of large differences 



compared to a pure quadratic function. So, we chose 
the Huber function [29] (see Fig. 1) 



if a (x) 



x 2 if \x\ < a 

2a\x\ — a 2 elsewhere 



which is convex and gives an acceptable modeling of 
the desired image properties. The a parameter tunes 
the trade-off between the quadratic and the linear part 
of the function. 



Fig. 1. Penalization functions ip: quadratic (lhs) and Huber 
(rhs). 



The criterion J~ Rcg is convex by construction and 
presents a unique global minimum: the optimization 
can be achieved by iterative gradient-like optimiza- 
tion techniques and we have implemented a pseudo- 
conjugate gradient procedure with a Polak-Ribieres 
correction method [30]. 



3.4 Optimization Stage 



The optimization process requires numerous evalua- 
tion of J/acg and its gradient hence numerous non- 
uniform DFT computations. In order to avoid these 
computations, J LS is rewritten, without changing the 
formulation of the problem. The new expression is 
founded on Toeplitz property of H^H and reads (see 
Appendix for details): 



L-l [ N-l 

:U/) = EN 2 -2K E /rUAv, 

1=0 I n,m=0 

N-l 

+ E C u ,v G u>v (6) 

u,v=l-N 

where C is the image correlation matrix, computable 
by FFT D and G are given by: 



D r 



G 1: 



1 L_1 

E s i e~ i2n ^ m+k y n ^ 



N 



1=0 

1 L - x 



N 2 



E 

1=0 



(7) 
(8) 



for n, m = 0, • • • , N — 1 and u, v = 1 — N, ■ ■ ■ ,N — 1 
and can be precomputed before the optimization stage. 

The 2N - 1 x 2N - 1 matrix G depends on the k- 
space trajectory only and can be computed once for 
all, given a trajectory. Moreover, it has a Hermitian 
symmetry, G^ = G, which allows to compute only 
one half of the matrix. The N x N matrix D depends 
on the A;-space trajectory and on the measured data. 
It can then be precomputed, but must be recomputed 
with each new data set. 

The new expression allows to reduce the computa- 
tional complexity of the optimization stage: instead of 
one DFT computation at each iteration, only one pre- 
computed DFT is required, the criterion and its gradi- 
ent can be computed from D and G by means of usual 
products and FFT. 

The gradient using a matrix formulation, is given then 
as (see also Appendix for details): 



dJJJ) 
df 



2f*G-2D. 



where * is a bidimentional convolution efficiently 
computed by FFT. 



4 Simulation and acquisition results 



In this section the proposed reconstruction method is 
compared to the gridding method on a mathematical 
model and a real phantom both acquired using a spiral 
sequence. 



4.1 Simulated model 



The simulated model is a 128 x 128 complex valued 
image and mimics two vessels on a variable back- 
ground. The magnitude image includes homogeneous 



regions and sharp transitions, while the phase im- 
age, related to the velocity image, corresponds to a 
parabolic and a blunt flow profile on a zero phase 
background (see Fig. 2). 




Fig. 2. Simulated phantom: magnitude image on the left- 
-hand side (lhs) and phase image on the right-hand side 
(rhs). We have selected two ROIs: ROI1 is the central 
square and ROI2 is the blunt flow area (upper right circle). 



For the direct problem, i.e simulating the acquired 
data, the exact model has been used without any ap- 
proximation, which allows to compute the value of the 
/c-space data along any sampling trajectory. A data set 
of 6 spiral arms of 512 samples each have been sim- 
ulated, thus the number of samples (6 x 512) was 5 
folds less than the number of pixels (128 x 128). The 
reduced number of samples and their very irregular 
density makes the reconstruction problem non invert- 
ible and thus allowed to test the quality of the regu- 
larized reconstruction in the case of sparse data. 

The hyperparameters, chosen empirically to obtain the 
best possible reconstruction, have been set to follow- 
ing values: Ai = 0.1, a\ = 20, A = 0.5 and a = 10 
and were then also used for the phantom reconstruc- 
tion. A 7 x 7 Kaiser-Bessel kernel, as introduced in 
[6], was used for the gridding reconstruction. 

It can be observed that the regularized reconstruction 
offers a better visual quality than the gridding (Fig. 
4) and that it is closer to the reference image. Sharp 
edges are maintained and enhanced while at the same 
time the noise level is smoothed throughout the im- 
age. This trade-off is achieved by the properties of 
the selected penalization function. The reconstruction 
presents less artifacts inside and outside the inner part 
of the image while the spatial resolution is preserved. 
These aliasing artifacts due to the undersampling are 
greatly reduced but their structure is more complex 
to analyze than for a Cartesian acquisition due to the 
characteristics of the spiral sampling trajectory [31]. 
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Fig. 3. Sensitivity to hyperparameters around the chosen 
values Ai = 0.1, a\ = 20, Ao = 0.5 and a® = 10: re- 
construction errors when one hyperparameter is varied at 
a time. Case for 6 spirals and 512 samples/spiral (noise 
free), selected values are indicated (as dots) on each curve. 
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Fig. 4. Reconstruction for noisy data (30 dB): 6 spirals and 
512 samples/spiral: re-gridding method (lhs) and proposed 
one (rhs). The top part shows the modulus images, middle 
part shows rows 50 and 75 and bottom one shows difference 
images with the reference. 



The examination of the A;- space of the reconstructed 
images (FFT of the reconstructed images), shown in 
Fig. 5, allows to compare the frequency content of the 
two reconstructed images versus the reference one. 
The proposed method restores a A;-space very close 
to the reference one, while the gridding reconstruc- 
tion still lets appear the underneath sampling trajec- 
tory. This shows that the a priori introduced by the 
regularization is more pertinent and helps to restore 
an image closer to the original object. 
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Fig. 5. From left to right : reference /c-space, gridding 
A;-space and reconstructed /c-space (6 spirals and 512 sam- 
ples/spiral). 

Figure 6 presents a quantitative validation of the 
method, varying the number of spirals, the number of 
samples per spiral and the SNR, using the following 
criteria. 

• The quadratic reconstruction error in ROI1 (see 
Fig. 2) which gives a measure of the distance be- 
tween the reconstruction and the reference. 

• The variance for the constant gray level region of 
ROE (see Fig. 2) [13]. 

These figures confirm the former qualitative results. 
The proposed method gives a quadratic error 5 to 300 
folds lower than the gridding, while the variance is 
improved 3 to 10 folds whatever the sampling or noise 
level. 

Figure 3 presents a quantitative evaluation of the hy- 
perparameter sensitivity computed as the variations of 
the squared reconstruction error in a defined region of 
interest (ROI1). We note that the selected values are 
very close to the ones that minimize the errors when 
only one hyperparameter is varied at a time. The in- 
tervals where these parameters can be chosen are rel- 
atively large: this ensures that the solution is robust 
with respect to the hyperparameter values. 

This results show that the quality of the image can 
be maintained while using acquisition sequences that 
sample a smaller number of data and then reduce the 



acquisition time, proportionally to the number of ac- 
quired spirals. 



4.2 Phantom acquisition 



The method was then tested on the GEMS test phan- 
tom with a 1.5 Tesla Signa systerrQ- The sampling 
trajectory consisted in 24 interleaved spirals each of 
2048 samples and a 16 cm FOV. 

Figure 7 presents the reconstructed magnitude images 
(256 x 256) for the gridding and the regularized meth- 
ods as well as a zoom in the comb like ROI for the gen- 
uine acquisition geometry (24 spirals). Fig. 8 presents 
the corresponding results when one spiral over two 
has been discarded, providing a gain of two in the ac- 
quisition time. 
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Fig. 7. Reconstruction with 24 spirals and 2048 sam- 
ples/spiral. From top to bottom : modulus image, ROI. On 
the left the gridding reconstruction and on the right the 
proposed method. 

For the genuine acquisition geometry (Fig. 7) the reg- 
ularized image is very close to the gridding recon- 
struction and it even shows a slight reduction of the 
noise level in the background. 

Undersampling strongly degrades the gridding recon- 
struction (Fig. 8): only a small central region remains 
free of all artifacts. As has been shown for simulated 



1 Acquisition are provided by MJ. Graves, University of 
Cambridge and Addenbrooke's Hospital, Cambridge, UK. 
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Fig. 6. First row: quadratic reconstruction error in ROI1 and second row : variance in ROI2. The gridding method results 
are plotted with a dotted line and the proposed method results with a solid line. From left to right: 10 spirals, variable 
number of samples (128, 256 et 512), no noise; variable number of spirals (4 to 10), 512 samples/spiral, no noise; 10 
spirals, 512 samples/spiral, variable noise level (20 to 50 dB). 
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Fig. 8. Reconstruction with 12 spirals and 2048 sam- 
ples/spiral. From top to bottom: modulus image, ROI. On 
the left the gridding reconstruction and on the right the 
proposed method. 

data, our method applied to actual measurements gives 
an image where the aliasing artifacts are strongly re- 
duced inside the object and where a very homoge- 
neous background is preserved. The two comb-like 
ROIs (Fig. 8) show more clearly the improvement pro- 
vided by the proposed method. The regularization also 
provides an image with well defined edges, illustrating 
that the chosen prior is well suited to achieve the com- 
promise between noise smoothing and contour preser- 



vation constraints. 

Characterization of aliasing artifacts can be ap- 
proached by studying the structure of the matrix G 
which can be interpreted as the point spread function 
of the imaging system: the observed image being the 
convolution of the true object with G. Fig. 9 shows 
this matrix for the two sampling schemes. The central 
white spot (resp. peak in the ID figures) introduces 
a blurring effect proportional to its diameter (resp. 
width), while the outer circles (resp. peaks) are re- 
sponsible for aliasing. The closer these circles to the 
center the more important the aliasing artifacts. The 
undersampling that shrinks these circles was partially 
inverted by the proposed method while it was kept 
unchanged by the gridding reconstruction. 



Beforehand computation of matrices D and G con- 
siderably speeds up the optimization procedure. How- 
ever, if G can be computed once for all for a given 
acquisition sequence, D must be computed for each 
data set. The computational complexity that arises in 
computing D is not a drawback for clinical use of the 
method: it takes 30 sec to compute matrix D (12 spi- 
rals, 2048 samples per spiral, image 256 x 256) using 
a C-Program on a PC computer with an AMD-Athlon 
2.1 GHz processor. 




Fig. 9. Matrix G for 24 spirals (lhs) and 12 spirals (rhs) 
(Log scale). Sharp peaks denoted by arrows cause aliasing. 

The optimization was performed using the computing 
environment Matlab in 3 minutes, and 50 iterations 
were needed to converge to an accurate solution. Each 
iteration requires one gradient and three criterion cal- 
culations. The calculations of the FFTs represent the 
main computational burden during the minimization: 
every iteration involves six 512 x 512 2D-FFTs for the 
criterion and two 768 x 768 2D-FFTs for the gradi- 
ent. This time could be considerably reduced by im- 
plementing the algorithm on a dedicated processor. 
Indeed, given the characteristics of the Texas Instru- 
ments TMS320C64x series DSP, all of the FFTs could, 
theoretically, be performed in about 18 sec, leading to 
an important decrease in the total optimization time. 

Moreover the computation of the criterion, the gra- 
dient and the matrix D are highly amenable to par- 
allelization, and with a sufficient number of process- 
ing elements, the reconstruction could be done even 
faster, which could allow the use of the method in a 
wide variety of clinical applications. 



5 Discussion / Conclusion 

The proposed method differs from more conventional 
ones insofar as it does not involve any regridding of 
the acquired data and accounts for edge preserving 
smoothing penalties. Utilization of only the acquired 
data and integration of smoothness and edge preserva- 
tion penalization in the reconstruction opens the way 
to strong improvement in MRI. 



From a computational stand point, the original for- 
mulation leads to the awkward situation of an opti- 
mization algorithm permanently shifting from Fourier 
to image domains requiring for numerous heavy non- 
uniform Fourier transform computations. Rewriting 
the criterion allowed to perform the whole opti- 
mization in the image domain providing the pre- 
computation of two matrices. The first one character- 
izes the geometry of acquisitions in A;-space and gives 
interpretation of aliasing structures; the second can be 
seen as a discrete Fourier transform of the acquired 
data. 

Moreover, alternatives exist to still improve the recon- 
struction efficiency of the method: substituting non- 
uniform FFT algorithms for the non-uniform Fourier 
transform in the pre-computations [32]; calculating a 
solution corresponding to a small ROI only; substi- 
tuting a Newton like [33] or a dual optimization [34] 
method to the conjugate gradient could dramatically 
reduce the computational cost and make the method 
available for clinical applications. 

Finally, the inverted model could be improved by in- 
tegrating an exponential term that takes into account 
the relaxation of the magnetic moments. A Laplace in- 
version framework should then be substituted for the 
present Fourier framework but the overall inversion 
procedure will remain valid. 
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Appendix 

The appendix gives detailed calculi for the new form 
of the fit to the data term <J LS and its gradient, required 
for efficient numerical optimization. 



A Criterion calculus 

We have: 

JUf) = E \si-yi\ 2 = E \si\ 2 +\ yi \ 2 -2K{ Siy n . 

1=0 1=0 

(A.l) 



where yi — hif is the noise free model output given by 
Eq (2). It is the sum of quadratic terms over the whole 
acquired data. The first term is simply the norm of the 
data, the second one is developed in subsection A.l 
and the last one in subsection A.2. 



A.l Term involving model output yi 



Expansion of \yi\ 2 , given model (2) yields: 
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A change of the summation variable : u = p — p' and 
v = q — q' gives: 

q=M(v) 
i N-l p=M(u) 
I |2 _ J_ V V f f* J2-K(k l x u+k l y v) 
~~ ]\[2 JP,q Jp-u,q-v K 



u,v=l—N p=m(u) 
q=m(v) 



where m(-) and M(-) are the index summation bounds: 
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after rearrangement of the summations. Let us state 

foru,v=l-N,...,N-l: 
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(A.2) 



which only depends upon the A;-space trajectory. We 
finally have: 

L-l N-l 

El^| 2 = E Cu,vG U;V . (A.3) 

1=0 v,u=l-N 



A.2 Term involving model output y\ and observed 
data si 



Using (2), the involved term writes: 
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and summation over I yields: 
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and 



M(w) = 



N-l 
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iV-l-|ii;| ifw<0 
We then introduce the image correlation matrix 



q=M(v) 
p=M(u) 

C u ,v = E fp>q /] 

p=m(u) 
g=m(i)) 
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which can be computed by FFT. So, |^| 2 simply 
writes: 
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The summation over / yields 
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p,q=0 1=0 



after rearrangement. We then introduce the DFT, for 

p,q = 0,...,N-l: 



Dp,q = ^Y. S l^ (kiP+klyq) ( A - 4 ) 
iV 1=0 

which depends upon observed data and A:-space tra- 
jectory. The current term then simply writes: 

5>!/r= E/; )9 A>,<r (A.5) 
z=o p,g=o 

Substitution of (A.5) and (A.3) in (A.l) yields the 
announced form of Eq. (6). 
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B Gradient of the Criterion 



The partial derivative of (A. 5) with respect to /„ 
clearly writes: 

Q L-l Q N-l 

s iUi = X] fp,q Dp,q = D n ^ m 

u Jl 



Of, 
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nm p t q=Q 



The partial derivative of (A. 3) with respect to f nm is 
more complicated. 
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Finally, we can write, using the expressions of the 
matrices D and G: 



d 7V " 1 



nm V>U=1 _ N 
L-l N-l 



- _ V V f -i2ir(k l x u+k l v) 
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where G* is the conjugate of G. 

The total gradient using a matrix formulation, is given 
then as: 



df 



2f*G-2D. 



where * is a bidimentional circular-convolution that 
can be efficiently computed by FFT 
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